Magnetostatics solution (nonlinear case)ΒΆ
[1]:
%%capture
%run current_density.ipynb
%run TEAM13_nonlinearity.ipynb
[2]:
##############################################################################
# Tree/Cotree gauging
##############################################################################
MESH = pde.mesh3.netgen(geoOCCmesh)
R = pde.tools.tree_cotree_gauge(MESH)
[3]:
R
[3]:
<27728x23597 sparse matrix of type '<class 'numpy.float64'>'
with 23597 stored elements in Compressed Sparse Column format>
[4]:
##############################################################################
# Assembly
##############################################################################
linear = '*coil,default'
nonlinear = 'r_steel,l_steel,mid_steel'
maxIter = 100
order = 1
fem_linear = pde.int.evaluate3(MESH, order = order, regions = linear).diagonal()
fem_nonlinear = pde.int.evaluate3(MESH, order = order, regions = nonlinear).diagonal()
D = pde.int.assemble3(MESH, order = order)
phix_Hcurl, phiy_Hcurl, phiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = order)
curlphix_Hcurl, curlphiy_Hcurl, curlphiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = order)
aJ = jx_L2 @ D @ phix_Hcurl.T +\
jy_L2 @ D @ phiy_Hcurl.T +\
jz_L2 @ D @ phiz_Hcurl.T
aJ = 1e7*aJ
def gss(A):
curl_Ax = curlphix_Hcurl.T@A; curl_Ay = curlphiy_Hcurl.T@A; curl_Az = curlphiz_Hcurl.T@A
Kxx = curlphix_Hcurl @ D @ sp.diags(fxx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fxx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphix_Hcurl.T
Kyy = curlphiy_Hcurl @ D @ sp.diags(fyy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fyy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiy_Hcurl.T
Kzz = curlphiz_Hcurl @ D @ sp.diags(fzz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fzz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiz_Hcurl.T
Kxy = curlphiy_Hcurl @ D @ sp.diags(fxy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fxy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphix_Hcurl.T
Kxz = curlphiz_Hcurl @ D @ sp.diags(fxz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fxz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphix_Hcurl.T
Kyx = curlphix_Hcurl @ D @ sp.diags(fyx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fyx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiy_Hcurl.T
Kyz = curlphiz_Hcurl @ D @ sp.diags(fyz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fyz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiy_Hcurl.T
Kzx = curlphix_Hcurl @ D @ sp.diags(fzx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fzx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiz_Hcurl.T
Kzy = curlphiy_Hcurl @ D @ sp.diags(fzy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fzy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear)@ curlphiz_Hcurl.T
return Kxx + Kyy + Kzz + Kxy + Kxz + Kyx + Kyz + Kzx + Kzy
def gs(A):
curl_Ax = curlphix_Hcurl.T@A; curl_Ay = curlphiy_Hcurl.T@A; curl_Az = curlphiz_Hcurl.T@A
return curlphix_Hcurl @ D @ (fx_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fx_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) +\
curlphiy_Hcurl @ D @ (fy_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fy_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) +\
curlphiz_Hcurl @ D @ (fz_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + fz_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) - aJ
def J(A):
curl_Ax = curlphix_Hcurl.T@A; curl_Ay = curlphiy_Hcurl.T@A; curl_Az = curlphiz_Hcurl.T@A
return np.ones(D.size)@ D @(f_linear(curl_Ax,curl_Ay,curl_Az)*fem_linear + f_nonlinear(curl_Ax,curl_Ay,curl_Az)*fem_nonlinear) -aJ@A
import time
A = np.zeros(curlphix_Hcurl.shape[0])
mu = 0.0001
eps_newton = 1e-5
factor_residual = 1/2
tm = time.monotonic()
for i in range(maxIter):
gssu = R.T @ gss(A) @ R
gsu = R.T @ gs(A)
tm = time.monotonic()
# wS = chol(gssu).solve_A(-gsu)
wS = sp.linalg.spsolve(gssu, -gsu)
print('Solving took ', time.monotonic()-tm)
w = R@wS
alpha = 1
# ResidualLineSearch
# for k in range(1000):
# if np.linalg.norm(gs(A+alpha*w),np.inf) <= np.linalg.norm(gs(A),np.inf): break
# else: alpha = alpha*factor_residual
# AmijoBacktracking
float_eps = 1e-12; #float_eps = np.finfo(float).eps
for kk in range(1000):
if J(A+alpha*w)-J(A) <= alpha*mu*(gsu@wS) + np.abs(J(A))*float_eps: break
else: alpha = alpha*factor_residual
A_old_i = A
A = A + alpha*w
print ("NEWTON: Iteration: %2d " %(i+1)+"||obj: %2e" %J(A)+"|| ||grad||: %2e" %np.linalg.norm(R.T @ gs(A),np.inf)+"||alpha: %2e" % (alpha))
# if ( np.linalg.norm(R.T @ gs(A),np.inf) < eps_newton):
# break
if (np.abs(J(A)-J(A_old_i)) < 1e-5):
break
elapsed = time.monotonic()-tm
print('Solving took ', elapsed, 'seconds')
Solving took 2.796000000089407
NEWTON: Iteration: 1 ||obj: -3.220290e+10|| ||grad||: 2.479687e+07||alpha: 6.250000e-02
Solving took 2.5620000001508743
NEWTON: Iteration: 2 ||obj: -6.401723e+10|| ||grad||: 5.435596e+07||alpha: 1.000000e+00
Solving took 2.672000000020489
NEWTON: Iteration: 3 ||obj: -6.639095e+10|| ||grad||: 5.773030e+07||alpha: 5.000000e-01
Solving took 2.6410000000614673
NEWTON: Iteration: 4 ||obj: -1.508296e+11|| ||grad||: 6.567249e+07||alpha: 1.000000e+00
Solving took 2.5780000002123415
NEWTON: Iteration: 5 ||obj: -1.530142e+11|| ||grad||: 4.925437e+07||alpha: 2.500000e-01
Solving took 2.547000000020489
NEWTON: Iteration: 6 ||obj: -1.858670e+11|| ||grad||: 3.979715e+07||alpha: 1.000000e+00
Solving took 2.6089999999385327
NEWTON: Iteration: 7 ||obj: -1.878539e+11|| ||grad||: 3.481956e+07||alpha: 1.250000e-01
Solving took 2.515000000130385
NEWTON: Iteration: 8 ||obj: -1.909770e+11|| ||grad||: 2.171232e+07||alpha: 2.500000e-01
Solving took 2.594000000040978
NEWTON: Iteration: 9 ||obj: -1.931830e+11|| ||grad||: 1.881516e+07||alpha: 5.000000e-01
Solving took 2.6089999999385327
NEWTON: Iteration: 10 ||obj: -1.933408e+11|| ||grad||: 2.100793e+07||alpha: 2.500000e-01
Solving took 2.5310000001918525
NEWTON: Iteration: 11 ||obj: -1.936156e+11|| ||grad||: 1.285842e+07||alpha: 1.000000e+00
Solving took 2.5
NEWTON: Iteration: 12 ||obj: -1.939955e+11|| ||grad||: 6.739643e+06||alpha: 2.500000e-01
Solving took 2.547000000020489
NEWTON: Iteration: 13 ||obj: -1.940739e+11|| ||grad||: 3.975133e+06||alpha: 5.000000e-01
Solving took 2.530999999959022
NEWTON: Iteration: 14 ||obj: -1.941111e+11|| ||grad||: 1.127736e+06||alpha: 1.000000e+00
Solving took 2.5
NEWTON: Iteration: 15 ||obj: -1.941138e+11|| ||grad||: 1.987664e+05||alpha: 1.000000e+00
Solving took 2.5
NEWTON: Iteration: 16 ||obj: -1.941138e+11|| ||grad||: 9.657019e+03||alpha: 1.000000e+00
Solving took 2.530999999959022
NEWTON: Iteration: 17 ||obj: -1.941138e+11|| ||grad||: 7.475946e+00||alpha: 1.000000e+00
Solving took 2.5469999997876585
NEWTON: Iteration: 18 ||obj: -1.941138e+11|| ||grad||: 1.351817e-02||alpha: 1.000000e+00
Solving took 2.547000000020489
NEWTON: Iteration: 19 ||obj: -1.941138e+11|| ||grad||: 2.625993e-05||alpha: 1.000000e+00
Solving took 2.5939999998081475 seconds
[5]:
curlphix_Hcurl_P0, curlphiy_Hcurl_P0, curlphiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = 0)
Bx = curlphix_Hcurl_P0.T @ A
By = curlphiy_Hcurl_P0.T @ A
Bz = curlphiz_Hcurl_P0.T @ A
[6]:
##############################################################################
# Storing to vtk
##############################################################################
grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_L2_Vector(grid,Bx,By,Bz,'grad_x')
pde.tools.vtklib.writeVTK(grid, 'magnetostatics_solution.vtu')
[7]:
import pyvista as pv
mesh = pv.read('magnetostatics_solution.vtu')
mesh2 = pv.read('current_density.vtu')
mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,4])
p = pv.Plotter()
p.add_mesh(threshed, style='surface', color = "w", opacity=0.2, label=None)
threshed.set_active_vectors("grad_x")
arrows = mesh.glyph(scale="grad_x", orient=True, tolerance=0.03, factor=18)
p.add_mesh(arrows, color="orange")
mesh2.set_active_vectors("J_L2")
arrows2 = mesh2.glyph(scale="J_L2", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows2, color="black")
p.camera_position = [(0, -500, 200),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend="html")
[ ]: